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Molecular dynamics simulations are used to study different definitions of contact at the atomic 
scale. The roles of temperature, adhesive interactions and atomic structure are studied for simple 
geometries. An elastic, crystalline substrate contacts a rigid, atomically flat surface or a spherical 
tip. The rigid surface is formed from a commensurate or incommensurate crystal or an amorphous 
solid. Spherical tips are made by bending crystalline planes or removing material outside a sphere. 
In continuum theory the fraction of atomically flat surfaces that is in contact rises sharply from 
zero to unity when a load is applied. This simple behavior is surprisingly difficult to reproduce 
with atomic scale definitions of contact. Due to thermal fluctuations, the number of atoms making 
contact at any instant rises linearly with load over a wide range of loads. Pressures comparable to the 
ideal hardness are needed to achieve full contact at typical temperatures. A simple harmonic mean- 
field theory provides a quantitative description of this behavior and explains why the instantaneous 
forces on atoms have a universal exponential form. Contact areas are also obtained by counting the 
number of atoms with a time-averaged repulsive force. For adhesive interactions, the resulting area 
is nearly independent of temperature and averaging interval, but usually rises from zero to unity 
over a range of pressures that is comparable to the ideal hardness. The only exception is the case 
of two identical commensurate surfaces. For nonadhesive surfaces, the mean pressure is repulsive if 
there is any contact during the averaging interval At. The associated area is very sensitive to At and 
grows monotonically. Similar complications are encountered in defining contact areas for spherical 
tips. Even for the adhesive case, the area based on time-averaged forces can not be described by 
continuum theory. 

PACS numbers: 46.55.+d, 62.20.Qp, 81.40.Pq 



I. INTRODUCTION 

The area of intimate contact between surfaces plays a 
central role in continuum models of friction and adhesion. 
Due to surface roughness, this real area of contact A reSi \ is 
typically much less than the apparent area of the surfaces 
Aq. Analytic and numerical work indicates that A rca x 
rises linearly with load for nonadhesive surfaces-whether 
they deform elastically [H-3 or plastically 0-0] • This 
linear relation, and the assumption of a constant shear 
stress, provides one of the most common explanations 
of the linear relation between friction and load in many 
experiments 

The advent of nanotechnology and the ability to mea- 
sure friction in contacts with molecular dimensions has 
led to great interest in extending continuum models 
to nanometer scales and identifying its limits [8-23]. 
One fundamental question is what contact means at the 
molecular scale. The mere presence of an interaction be- 
tween surfaces is not enough, since attractive van der 
Waals interactions extend to arbitrary scales. One com- 
mon approach has been to associate the onset of direct 
repulsion between atoms with contact [l2l - [l5l . [l8l H3-H1] ■ 
For most systems this coincides with common definitions 
of atomic diameters [27| , but the notion of contact is still 
complicated by thermal fluctuations and the finite range 
of interactions. These are not commonly included in con- 
tinuum models and recent work reveals that they greatly 
complicate the extension of continuum views of contact 
to atomic scales. 



Early experiments found that friction measurements 
could be fit by applying continuum theory to nanometer 
scale contacts and assuming a constant shear stress at 
the interface [sT-flll [28[ . Simulations with tips of the same 
size showed that while friction forces could be fit in this 
way, the fit parameters were often significantly different 
than the independently determined material properties 
[l2l [l3| . The actual area of contact was often a factor 
of two larger than the continuum prediction, the contact 
stiffness was up to an order of magnitude smaller, and 
the friction force varied by two orders of magnitude with 
changes in atomic geometry by much less than an atomic 
diameter. 

Knippenberg et al. fl8j examined contact between 
atoms in sliding contacts between nanometer scale tips 
and a substrate covered with a surfactant. They found 
that the area of contact was broadened due to the com- 
pliance of the surfactant layer, but that most of the 
atoms contributed little to the normal and lateral (fric- 
tion) forces. Roughly 90% of the force was carried by 
only 10% of the atoms, calling into question the binary 
nature of contact. Cheng et al. [26| considered the effect 
of single monolayers of short chains. They also found 
very large variations in the force on individual atoms, 
both spatially and temporally. One consequence was that 
measures that included the magnitude of forces, such as 
moments of the pressure distribution, gave very different 
contact radii than simply counting the atoms feeling any 
nonzero repulsion. 

Mo et al. have considered contact between a bare 
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diamond substrate and an amorphous surface that is 
nominally spherical [3 UH . They found that the num- 
ber of contacting atoms grew linearly with the applied 
load. They concluded that their surfaces were sufficiently 
rough that the continuum results predicting area propor- 
tional to load could be applied. This idea was tested by 
comparing the distribution P(f) of the magnitude of lo- 
cal forces / to Persson's continuum theory |3| . 

Yang et al. [HI, [23| have tested Persson's theory by 
simulating contact between a flat substrate and a sur- 
face with roughness on all length scales. They found 
many more atoms with low forces than predicted. The 
area of contact was not obtained by counting atoms, but 
by fitting the distribution P(f) at larger forces. The 
result was then in reasonable agreement with Persson's 
predictions. Luan et al. [2|| [29|, [3(| examined two dimen- 
sional models with much larger linear dimensions. They 
also found many atoms had low forces and that P(f) 
decreased monotonically with /. The contact area ob- 
tained by counting atoms was extremely sensitive to the 
detailed atomic structure. One might not expect that 
continuum theory could capture these local details, but 
it did correctly capture the long-range stress correlations 
in the simulations [HI, [3lTl33l | . 

In this paper we explore the effect of thermal fluctu- 
ations on different definitions of contact. The simplest 
possible case, contact between atomically flat surfaces, is 
considered first. Even identical crystals with atoms di- 
rectly over each other exhibit unexpected behavior. The 
number of atoms exhibiting repulsion at any instant rises 
linearly with load, as might be expected for rough sur- 
faces in continuum theory. This linear scaling is ob- 
served for temperatures from 10~ 4 to 1/4 of the melt- 
ing temperature, T m . The contact area determined from 
time-averaged forces shows full contact for this geome- 
try, but not for incommensurate crystals with different 
lattice constants or surfaces cut from amorphous solids. 

A simple harmonic mean-field model for thermal fluc- 
tuations of individual surface atoms provides a quanti- 
tative description of the above results. It also explains 
several other new observations in our simulations. These 
include a nearly universal exponential distribution of the 
magnitude of instantaneous repulsive forces, and a cor- 
responding universal distribution of the fraction of load 
carried by the atoms with the largest force. There is 
also a direct correlation between the fraction of time in 
contact and the time averaged force on an atom that is 
independent of surface structure. The force required for 
atoms to remain in contact more than half the time is 
surprisingly large at typical temperatures, with the cor- 
responding pressure comparable to the ideal hardness. 

The prototypical case of contact between a spherical 
tip and flat substrate is considered next. Different atomic 
geometries of the tip, including commensurate, incom- 
mensurate and amorphous are considered. The same har- 
monic model used for flat surfaces describes the time av- 
erage force and fraction of time in contact for tip atoms. 
For nonadhesive tips, the contact area obtained by count- 



ing the number of atoms with a repulsive time-averaged 
force depends on the observation time and grows mono- 
tonically. Results for adhesive surfaces are less sensitive 
to time interval and temperature, but differ substantially 
from continuum theory. 

The paper is organized as follows. In Sec. II we de- 
scribe our simulation methods. Then MD results for the 
contact area between flat surfaces are described in Sec. Ill 
and a model for the distribution of contact forces is devel- 
oped and tested. Section IV examines contact between a 
sphere and flat and discussions and conclusions are pre- 
sented in Sec. V. 



II. SIMULATION METHODS 

In continuum theory, contact between two rough elas- 
tic solids can be mapped to that between a flat elastic 
substrate and a rigid rough upper solid @. We consider 
the latter case, and take the upper solid to be atom- 
ically flat or a spherical tip. Previous studies at zero 
temperature indicate that the continuum mapp ing re- 
mains approximately valid at atomic scales [12j . but it 
does not take into account thermal effects. At any finite 
temperature, annealed roughness from thermal fluctua- 
tions will be superimposed on top of any quenched struc- 
tural roughness on the surfaces. While the mapping may 
not accurately capture the quantitative effect of thermal 
fluctuations on both surfaces, using a rigid upper sur- 
face and a flat elastic substrate minimizes the parameter 
space to be explored. Preliminary results for two elastic 
solids show the same effects and trends are mentioned 
where relevant. 

Atoms in the upper solid interact with those in the 
substrate via a truncated and shifted Lcnnard- Jones (L J) 
potential [34j : 

V(r) = 4e[(a/r) 12 - (a/rf - (a/r c ) 12 + (a/r c f] , (1) 

where r c is the cutoff length, e is the binding energy, and 
a is the atomic diameter. We express all physical quan- 
tities in terms of LJ units based on e, <r, and the mass 
rn of the substrate atoms. For example, the unit of time 
is r = y/ ma 2 / e. Where possible we plot dimensionless 
quantities in order to facilitate comparison with experi- 
ments and theoretical models. 

To model a nonadhesive contact, we take r c = 2 1 / 6 cr 
so that the LJ potential is purely repulsive. For adhesive 
contact, we use r c = 2.2<r and reduce the binding energy 
to 0.5e to make adhesion between surfaces weaker than 
the internal cohesion in the substrate. This prevents the 
creation of cracks in the substrate under negative load 
(separating surfaces). 

The substrate is a fee crystal with a (001) surface. To 
make the substrate as elastic as possible, nearest neigh- 
bor atoms in the substrate interact through a harmonic 
potential: 

V lJ {r)= 1 -k{r~df , (2) 
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where k is a spring constant and d is the equilibrium 
spacing between nearest neighbors. We take d — 2 1 / 6 <7 
to match the position of the minimum in the LJ po- 
tential in Eq. [TJ The spring constant k = 57e/cr 2 also 
matches the second derivative of the LJ potential at its 
minimum. Rather than listing the specific bonded pairs, 
spring forces are calculated for all atoms that are within 
a distance of 1.3er. This limits interactions to nearest 
neighbors, but does allow plastic deformation to occur 
when forces are extremely high. We limit our simula- 
tions to loads below this point. 

The equilibrium density of the substrate at k^T/e = 
is p = 1.0m /a 3 . The initial state is a crystal of this 
density. Periodic boundary conditions with period l x — 
£ y = 190.49(T are applied along the surface. The crystal is 
l z = 189.69cr high along the z direction and the bottom 
layer is held fixed to mimic the support that balances the 
external load. This depth is large enough to approximate 
contact of a sphere and a semi- infinite substrate [l3L f35j j . 
The substrate has the same fee structure and periodic 
boundary conditions in all simulations. As the temper- 
ature T changes from 10 _4 e/fcB (called the low T case) 
to T = 0.175c/A:b (called the high T case), the height of 
the substrate shrinks by about 0.3% because of a slight 
anharmonicity induced by changes in orientations of the 
springs. This slight breaking of cubic symmetry does not 
affect our results [26| . 

We studied three atomic geometries for rigid flat sur- 
faces. Two contain a single (001) plane of a fee crystal 
with a lattice constant d! . One is a commensurate surface 
with d' jd— 1. This surface is in perfect alignment with 
the substrate, with its atoms directly above the equilib- 
rium positions of atoms in the top layer of the substrate. 
The second is an incommensurate surface with a lattice 
constant d'/d = 1.12342. The third surface has an amor- 
phous structure. It is a thin sheet with thickness ~ 0.8cr 
cut from an amorphous solid with density 1.0m/er 3 . The 
thickness was chosen to produce roughly the same num- 
ber of surface atoms as the commensurate case. 

The spherical tips have radius R = 100(7 and four types 
of atomic geometries studied previously [T3, 0, 26]. The 
commensurate tip is made by bending a (001) plane of 
a fee crystal with lattice constant d! — d into a spheri- 
cal shape. The incommensurate tip is bent from a (111) 
plane of a fee crystal with d'/d = 1.12342. The stepped 
tip is obtained by carving a spherical shell out of a com- 
mensurate fee crystal with d! — d. The detailed struc- 
ture of the stepped tip depends on realization and in our 
case the bottom layer is a (001) face with 104 atoms. 
The amorphous tip is a spherical shell carved out of the 
same amorphous solid from which the flat amorphous 
surface is cut. Tips used in AFM experiments are likely 
to be closest in structure to the amorphous or stepped 
tips. The spherical tips and flat upper surfaces all have 
a finite thickness. This is irrelevant for simulations of 
nonadhesive contacts. For adhesive contacts, the finite 
thickness merely reduces the effective adhesion between 
two surfaces and does not affect the trends and conclu- 



sions presented below. 

The simulations are performed using the Large- 
scale Atomic/Molecular Massively Parallel Simulator 
(LAMMPS) developed at Sandia National Laboratories. 
This classical MD code utilizes spatial decomposition to 
parallelize the computations. Forces are calculated with 
the help of neighbor lists. A velocity- Verlet algorithm 
with a time step dt — 0.005t is used to integrate the 
equations of motion. The substrate is held at a fixed tem- 
perature T using a Langevin thermostat. The Langevin 
damping rate T is typically O.lr -1 , but results reported 
here are not sensitive to T up to at least T ~ O.Ht^ 1 . 

To illustrate the role of temperature, we report results 
for T = 10~ 4 e/fcB and 0.175e//cB- In some cases, we also 
report results for T — 0.07e/fce- Note that the melting 
temperature of the substrates would be T m ~ 0.7e/fcs if 
LJ interactions with the same length and stiffness were 
used. Thus the three temperatures above correspond to 
roughly j^T m , ^T m , and \T m . MD simulations ig- 
nore quantum effects that reduce thermal fluctuations. 
Simulations below 5 to 10% of T m are usually not rep- 
resentative of the behavior of real materials, but may be 
useful for illustrating trends. 

In most simulations, a constant normal load L is ap- 
plied to the top surface and the system is allowed to 
equilibrate before data are collected. The time for stress 
equilibration is a small multiple of the time for sound to 
propagate across the system ~ 20r. We equilibrated the 
system for at least 250r after each change in load, and 
data were typically averaged over a subsequent 500t. 

For contact between nominally flat surfaces, the nat- 
ural dimensionless measure of load is L/AqE*, where 
Aq = t x t v is the nominal area of the contacting surfaces, 
L/Ao is the mean pressure in the contact, and E* is the 
effective modulus. The latter is related to Young's mod- 
ulus E and the Poisson ratio v by E* = 25/(1 — v 2 ) and 
E* = 63e/cr 3 for our substrate (26|. The ratio L/A E* 
is approximately equal to the average strain along the 
z direction far from the substrate. For contact with a 
tip of radius 2?, the substrate dimensions are effectively 
infinite and the natural dimensionless measure of load is 
L/R 2 E*. 



III. CONTACT BETWEEN NOMINALLY FLAT 
SURFACES 

A. Nonadhesive Contact Area 

Molecular dynamics simulations provide complete in- 
formation about forces and positions, but it is not obvi- 
ous how this information should be used to identify the 
real area of contact. For nonadhesive contact, it is nat- 
ural to say an atom is in contact when it feels a force 
from the opposing surface. This can be generalized to 
adhesive contact by saying atoms separated by less than 
some distance are in contact. In the following we use 
the same distance as for nonadhesive contact, so atoms 
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are counted as contacting when the force between them 
is repulsive. This definition has also been used in dis- 
cussing AFM experiments [24| , and previous simulations 
[3 E, E® [M US HI]. One can define an area of 
real contact, A c , by multiplying the number of contact- 
ing atoms N c by the area per atom A a , i.e., A c = N c A a . 
Thermal fluctuations introduce ambiguities in this defi- 
nition because atoms fluctuate in and out of contact. 

We first consider the case of nominally flat non- 
adhesive surfaces. From the continuum perspective one 
would expect that these surfaces should be in complete 
contact, A c = Aq, at any positive load. However surfaces 
constructed of discrete atoms can never be perfectly flat. 
Even when all the rigid atoms are at the same height, as 
for commensurate and incommensurate cases, the height 
at which substrate atoms feel a given force varies with 
lateral position. This corrugation depends on the rela- 
tive size of substrate atoms and the normal force, but the 
peak-to-peak height change is of order cr/3 for the cases 
considered here [H, H3| ■ The change in height of atoms 
on the amorphous surface leads to additional roughness 
of the same order. Because the potential changes very 
rapidly with separation on these scales, the atomic scale 
corrugation of any surface can lead to large changes in 
contact properties even in the limit of zero temperature 
considered in Refs. [HIE]. 

Thermal fluctuations lead to additional time varying 
or annealed roughness on the surface. An estimate of 
the magnitude of local fluctuations can be obtained from 
the root mean squared (rms) normal displacement of a 
surface atom if its eight neighbors are held fixed: 

Sz lms « k B T/k cS , (3) 

where k c g — 2k. As ksT/e increases from 10~ 4 to 0.175, 
this increases from about 10 _3 cr to 0.04<t. These values 
are similar to the measured height fluctuations discussed 
further below. One can estimate the period of norm al 
vibrations T v n, in the same way: T v n, ss 27ry m/k e s « 
0.6t. 

While Szrms is always smaller than the quenched sur- 
face corrugations described above, thermal fluctuations 
lead to time variation that significantly complicates the 
definition of contact. In particular, the number and iden- 
tity of atoms that exert a force on the opposing surface 
varies with time. Some authors have associated A rea i 
with the mean number of contacting atoms at any instant 
in time [13, E, E]- However, continuum theories gen- 
erally consider time-averaged forces and displacements. 
From this perspective it may be more natural to find the 
average force on surface atoms over a long time inter- 
val and identify all atoms with an average repulsion as 
in contact. One can interpolate between these limiting 
cases by defining A c (At), the mean area associated with 
atoms that exert a repulsive force over an averaging in- 
terval At, where At varies from a single time step to the 
length of the run (2(|. Note that this quantity does not 
measure any time evolution or aging of the contact, but 
rather the mean number of contacting atoms over a fixed 



time interval averaged over all starting times in a steady 
state simulation. 

Fig.[T]shows A c (At)/Ao for commensurate, incommen- 
surate, and amorphous surfaces at k^T/e = 0.175 and 
L/AqE* = 5.5 x 10~ 4 . In all cases, only a few percent 
of the substrate atoms contact the rigid surface at any 
instant of time. This percentage is nearly constant for At 
much shorter than T v n, because atoms have not had time 
to move [26( | . Over longer averaging intervals, A c rises 
and appears to saturate. For the commensurate case, 
where all substrate atoms are directly under rigid atoms, 
A c saturates at full contact in a time of order 10r. For 
the other surfaces there is a slow, roughly logarithmic, 
rise over three decades in time before A c appears to sat- 
urate. The final value from the average over the entire 
simulation (5 x 10 4 r) is over 90% for the incommensurate 
surface, and only 30% for the amorphous case. 

The difference in the asymptotic values of A c reflects 
the atomic structure of the surfaces. The commensurate 
case is the closest to ideally flat since all substrate atoms 
are directly below rigid atoms and have the same mean 
spacing. As a result, all atoms are likely to contact within 
a time of order the vibrational time. The above estimate 
oiTyii, ~ 0.6r was for short wavelength modes. While A c 
rises rapidly as the time interval becomes of this order, 
complete saturation does not occur until times of order 
the period of the slowest long wavelength modes of the 
system ~ £ z /c s ~ 25r, where c s is the speed of acoustic 
waves. 

For the incommensurate case, each substrate atom has 
a slightly different environment. Some are directly below 
rigid wall atoms, and some are centered in between wall 
atoms. As noted above, this changes the height needed 
to produce a given force by ~ 0.3c This is larger than 
the estimated thermal height fluctuations, Sz lms ~ 0.04er. 
However, as At increases, more and more unlikely config- 
urations will be sampled and the total number of atoms 
that come into contact will grow. In principle, all atoms 
should contact given an arbitrarily long time, but over 
the times accessible to simulations the fractional contact 
area grows roughly logarithmically and then appears to 
saturate. The situation for amorphous walls is similar, 
but the rougher surfaces lead to a smaller long time value. 

Results for other temperatures and loads show qual- 
itatively similar behavior, but with different values of 
A c at early and late times. To illustrate trends with T 
and L, we will compare the mean results for atoms con- 
tacting over three time intervals: instantaneous ^4 C (0), a 
time interval At = 0.5t that is of order T v n,, and a time 
comparable to the entire simulation At = 500r (~ Ins). 
Fig. [2] shows results for commensurate, incommensurate 
and amorphous walls at two very different temperatures: 
ksT/e = 0.175 (open symbols) and 10~ 4 (closed sym- 
bols). The data are plotted on a log-log scale to capture 
the wide range of values. Note that averaging over even a 
single vibrational period of the atoms has a dramatic ef- 
fect on A c /Aq, leading to an order of magnitude increase 
in many cases. 
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FIG. 1: (Color online) Fractional contact area A c /Ao vs the 
time interval At on a log scale for three types of upper sur- 
faces: commensurate (solid line); incommensurate (dashed 
line); amorphous (dotted line). Here feT/e = 0.175 and 
L/A Q E* = 5.5 x 10" 4 . 



At small loads, all curves are nearly linear, implying 
that A c grows as a power of load. The lines show fits to 
a model developed in Section UlI CI that predicts a nearly 
linear relation between area and load. Best fits to the 
data for A c (0) at high T are consistent with this, giving 
exponents a little greater than 0.9 over several decades. 
Quenched geometrical disorder becomes more important 
at low temperatures, leading to changes in the apparent 
power law exponent for the amorphous and incommen- 
surate cases. 

For the commensurate surface, the area rises roughly 
linearly with load for all temperatures and time inter- 
vals. However the curves shift up with increasing At 
since more atoms have time to contact. There is a corre- 
sponding downward shift in the load needed to reach full 
contact. Over a time of 500r (~ Ins), nearly all atoms 
contact even at the lowest load and highest temperature. 
This full contact is consistent with one's expectations for 
perfectly fiat surfaces in continuum theory. 

Lowering the temperature increases the number of con- 
tacting atoms and lowers the load needed to achieve full 
contact. The rigid surface is held up by interactions 
with the highest substrate atoms, because the repulsion 
changes rapidly with separation. The roughness from 
thermal fluctuations grows with temperature as does the 
force from impacts at the thermal velocity. These effects 
lead to a decrease in A c with increasing T at a fixed load 
that is captured by the model discussed below. 

The incommensurate and amorphous surfaces show 
similar trends at high temperatures, except that A c /Aq 
saturates below unity on the time scale of our simu- 
lations. The saturation value is always lower for the 
rougher amorphous surface, but grows more rapidly with 
load. At low T, the variation with time interval is much 
weaker, with A c /Aq changing less than a factor of 2 with 



At at all loads for the amorphous case. The reason is 
that the quenched structural variations due to the incom- 
mensurate spacing or random atomic heights are large 
enough that some atoms are very unlikely to be brought 
into contact by thermal fluctuations at low T. Those sub- 
strate atoms that are closest to rigid atoms make contact 
quickly and those farther away never contact. 

It is interesting to compare the above results to a sim- 
ple model motivated by contact between a solid wall and 
ideal gas with number density n. When averaged over a 
long time, there is a uniform pressure p — nk B T on the 
wall from collisions with the gas. From the continuum 
perspective this implies that the entire solid surface is 
in contact with the gas. However at any instant, a very 
small fraction of solid atoms feel a very large collision 
force. The force is typically estimated from the change 
in momentum, using a typical thermal velocity normal 
to the wall v t = ^Jk B T /m. The momentum change in an 
elastic collision is 2mv t . Taking the vibrational time as 
an estimate of the collision time, the average force exerted 
on contacting atoms is ~ 2rn,v t /T v n 1 = vfeefffeT '/it. The 
fraction of area in contact will be the equilibrium pres- 
sure times the area per atom A a divided by this force. 
For typical values, the fractional contact area from this 
definition is extremely small even though the entire sur- 
face contacts the gas over long time intervals. 

For nonadhesive solid-solid contacts with short range 
interactions and high temperatures, the fraction of time 
in contact is also small. We can apply the same picture 
with the load replacing the ideal gas pressure. Then 

A c _ LA a 7T L 
A Q ~ A ^k cS k B T ~ CA A E* 

with the constant of proportionality 

irA a E* la 2 k eS ... 
o-fceff y k B T 

The first factor is of order unity since the same springs 
determine both E* and k e g. For k B T/e = 0.175, ca ~ 56 
which is comparable to the instantaneous values observed 
for all surfaces. The slightly slower than linear increase in 
A c /Aq with load is captured by the more detailed model 
developed in Sec. IIIC. 

Decreasing k B T/e from 0.175 to 10~ 4 should increase 
ca by a factor of 42. This is consistent with the data for 
commensurate surfaces. The increase is slightly smaller 
for incommensurate surfaces and even smaller for amor- 
phous surfaces. In contrast to the case of ideal gases, the 
quenched disorder on these surfaces keeps some atoms 
from touching even at large times and the number of 
these atoms grows as T decreases. As noted above, MD 
simulations ignore quantum effects that usually become 
important below 5 to 10% of the melting temperature. 
Thus thermal fluctuations should be even smaller at this 
low T, or equivalently it may be representative of the 
behavior at higher temperatures when quantum effects 
are included. Over the range of T where MD simulations 
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are accurate, the behavior is close to that for the high 
temperature in Fig. [2J 

To place the range of loads in perspective, the mean 
pressures in the contact, L/Aq, correspond to 0.5MPa 
(a few atmospheres) to 2GPa for a typical metal with 
E* ~ 200GPa. The maximum normal pressure that can 
be supported before plastic deformation is the hardness, 
H, and for bulk metals H/E* is typically 1(T 4 to 10" 3 . 
While significantly larger dimensionless hardnesses are 
observed in nanocrystals and in amorphous materials, 
Fig. [5] spans the range of dimensionless loads that are 
likely to be found in real materials. 
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FIG. 2: (Color online) Contact area A c vs load L for contact 
between a nonadhesive elastic substrate and a rigid flat upper 
surface with different geometries: (a) commensurate; (b) in- 
commensurate; (c) amorphous. Open and filled symbols are 
for T = 0.175e/fce and 10~ 4 e/fcB, respectively. The contact 
area is measured by counting the number of atoms in the top 
layer of the substrate that interact with the opposite surface 
at any instant At = (O) or during time intervals At = 0.5r 
(□) or 500r (A). The lines in the panel come from the model 
developed in Sec. III. C, with solid lines and dashed lines for 
fcsT/e = 0.175 and 10 -4 , respectively. 



B. Adhesive Contact Area 

The notion of contact is more complicated when long- 
range adhesive interactions are present. In principle, sur- 
faces feel a van der Waals interaction at arbitrarily large 
separations, although it will be arbitrarily weak. Some 
threshold for interactions must be introduced to create 
a sharp criterion for the onset of contact. A common 
and reasonable choice is the onset of repulsive interac- 
tions [HI. This corresponds to a separation distance at 
the minimum of the interatomic potential, which is com- 
monly used to define the atomic diameter. It is also the 
same as the separation used to define contact in the non- 
adhesive case for the LJ interactions used here. Small 
changes in this criterion have little effect on the trends 
described below. 

Fig. [3] shows the contact area A c (At) for adhesive in- 
teractions between commensurate, incommensurate and 
amorphous surfaces, at different time intervals and tem- 
peratures. A very strong adhesive potential correspond- 
ing to half the internal cohesive potential is used to max- 
imize the contrast with the nonadhesive case. This pro- 
duces very large tensile strains of a few percent that could 
produce yield if defects like dislocations were present. 
The range of interactions is extended to r c — 2.2ct, and 
further extension had little effect on the results. 

In all cases, adhesive interactions bind the surfaces in 
a local free energy minimum. The surfaces remain locked 
in this minimum until the load exceeds a negative thresh- 
old —Lp, where L p is often called the pulloff force in the 
context of tip-substrate interactions. The pulloff force is 
only unique in the limit of zero temperature, where it 
represents the load where the energy minimum becomes 
linearly unstable. At any finite temperature, thermal ac- 
tivation will eventually lead to pulloff at any constant 
negative load [Hj]. The pulloff force observed in simu- 
lations decreases with increasing observation time and 
temperature, and increases with surface area. 

As in simulations for tips [l3[ , the pulloff force is largest 
for the commensurate case where all substrate atoms are 
equally spaced from rigid atoms and can exert the max- 
imum adhesive force. The pulloff force is lower for the 
incommensurate case where all local environments are 
sampled, and even lower for the amorphous case where 
additional roughness reduces the number of atoms that 
contribute to the adhesion. This effect is entirely geo- 
metric since the interaction potentials are the same in all 
cases. 

The behavior of A c at high temperatures is similar to 
that for nonadhesive surfaces, although shifted to neg- 
ative loads. The fractional instantaneous contact area 
A c (0) rises roughly linearly at low fractional contact ar- 
eas for the commensurate case, with increasing negative 
curvature as area and load increase. The curvature is 
more pronounced for the other surfaces, but note that 
A c /Aq is bounded by unity and must saturate at large 
loads. The linear scale of Fig. [3] accentuates the large 
values of A c /Aq, where there was also more curvature for 
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FIG. 3: (Color online) Contact area A c vs load L for contacts 
between an adhesive elastic substrate and a rigid flat upper 
surface with different geometries: (a) commensurate; (b) in- 
commensurate; (c) amorphous. Open and filled symbols are 
for ksT/e = 0.175 and 10 -4 , respectively. The contact area is 
measured by counting the number of atoms in the top layer of 
the substrate that feel a repulsion from the opposing surface 
at any instant At = (O) or a t any point during an interval 
of At = 0.5r (□) or 500r (A). 



the nonadhesive case. Thermal fluctuations in the pullof 
force prevent us from accessing much smaller values of 
A c /Aq in constant load simulations, but the asymptotic 
behavior is discussed further in Sec. IIIC. 

For larger values of At, substrate atoms are counted as 
in contact if they feel a repulsive force from a rigid atom 
at any instant during the time interval. This is consistent 
with the definition for nonadhesive surfaces and leads to 
a monotonic rise in A c with At. For the commensurate 
case, atoms are very likely to make contact even within 
0.5r and full contact was obtained over the entire sim- 
ulation. The incommensurate surface also reached full 
contact for most loads at the longest time interval, while 
the greater roughness on the amorphous surface led to 
saturation at lower A c /Aq. 

At low temperatures, thermal fluctuations are smaller 



and the time interval is less important. For the commen- 
surate surface there is a sharp transition between zero 
and full contact at L/A E* ~ —0.01. The transition 
moves to negative loads as At increases because there is 
time for more distant atoms to contact. As for nonad- 
hesive surfaces, At has much less effect for incommen- 
surate and amorphous surfaces, because thermal fluctu- 
ations are smaller than the quenched variation in atomic 
separation. 

The above approach is only one way of generalizing 
the definition of contact used for nonadhesive surfaces. 
One could also say that atoms contact only when the 
time-averaged force over At is repulsive. For nonadhesive 
surfaces the force is always repulsive so the time-averaged 
force is repulsive if contact is made at any instant during 
St. For adhesive surfaces, intervals of repulsion can be 
countered by attractive interludes. From the continuum 
perspective, it may be most natural to define contact 
based on the time-averaged force. 

Fig. |4] shows the fraction of area that feels an average 
repulsion over time intervals of 0.5t and 500r. The in- 
stantaneous results are the same as in Fig. [31 Averaging 
reduces the thermal fluctuations about the mean force. 
For the amorphous and incommensurate surfaces, even 
averaging over At — 0.5r is sufficient to remove most 
of the fluctuations. The quenched disorder then domi- 
nates the variation with load and the results are nearly 
independent of both T and At. 

There is no quenched disorder for the commensurate 
case, so all atoms have the same mean force. For a suf- 
ficiently long time average, there is a sharp transition 
from no contact to full contact at L = 0. Shorter time 
averages give a Gaussian distribution of mean forces and 
A c /Aq rises like an error function. Both the Gaussian 
and the rise in A c /Ao sharpen as At grows. Note that 
the width of this rise for 0.5r and a commensurate sur- 
face, is smaller than the width of the rise for the other 
surfaces and large At. This is consistent with the lack of 
variation with At for these surfaces. 

The above results indicate that defining contact area 
based on time average forces provides a less ambiguous 
measure of contact for adhesive systems. Except in the 
extreme case of aligned, commensurate surfaces, the re- 
sults are insensitive to At. In contrast, measures based 
on instantaneous contact show an increasing drift with 
At and greater variation with T. The lack of negative 
forces in the nonadhesive case means that time average 
forces are sensitive to rare events. This difficulty can 
be overcome by averaging the position rather than the 
force. Unlike the force, the position is not positive def- 
inite. However, since the fraction of atoms in contact is 
almost always much less than a half for the loads con- 
sidered in the figures above, criteria based on the mean 
position almost always give zero contact area even though 
the local pressures are very high. 

Note that the association of repulsion with contact is 
also somewhat arbitrary. One could argue that surfaces 
that are bound together in a free energy minimum must 
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FIG. 4: (Color online) Contact area A c vs load L for con- 
tacts between an adhesive elastic substrate and a rigid flat 
upper surface with different geometries: (a) commensurate; 
(b) incommensurate; (c) amorphous. The contact area is de- 
termined from the mean number of atoms in the top layer of 
the substrate that feel a time-averaged repulsive force from 
the opposing surface over time intervals of 0.5r (O an d □) 
or 10 3 r (A and x) at k B T/e = 0.175 (O and A) or 10~ 4 (□ 
and x). 



be in contact. In this case A c /Aq would be significant at 
any L > —L p , while Figures [3] and 0] show zero contact 
area over much of this range. This could be addressed by 
associating the distance for contact with a larger spac- 
ing. For example, the commensurate surface becomes 
mechanically unstable when the atoms are separated by 
the distance r ~ 1.244a that maximizes the attractive 
LJ force. 



C. Force Distributions 

The results described above show that thermal fluctu- 
ations lead to large variations in the identity of atoms 
making contact and the force exerted by these atoms. 
Previous work has also shown that load may be very un- 



evenly distributed between atoms, with a very small frac- 
tion of the atoms carrying the bulk of the load [H, [26[ . In 
this section we examine the distribution of forces in time 
and among atoms. The distributions are surprisingly uni- 
versal and can be understood from a simple model that 
treats each atom independently. 

The distribution of forces at any instant can be charac- 
terized by the number of atoms with forces greater than 
some threshold /{, N(f t ), and the total force from these 
atoms F(f t ). Since we associate contact with repulsive 
(positive) forces, the discussion will focus on f t > 0, but 
it is readily extended to negative values. As f t decreases 
to zero, N(ft) and F(f t ) approach the total number of 
contacting atoms and the total repulsive force at that 
instant, respectively. The fraction of contacting atoms 
with / > f t , N(f t )/N(0), and the fraction of the re- 
pulsive force they carry, F(f t )/F(0), both vanish as f t 
increases. 

As shown in Figure [5] we find that the fraction of 
atoms carrying a given fraction of the total repulsive 
force is surprisingly universal. Results for all surfaces, 
with and without adhesion, have nearly the same form 
at k-aT/e = 0.175. Similar results are obtained at tem- 
peratures up to two orders of magnitude lower. Only 
at extremely low temperatures, such as 10 _4 e//cB, does 
the amount of quenched disorder become more important 
than thermal fluctuations. Even in this limit the amor- 
phous results remain similar to the high temperature re- 
sults. Results for the perfectly ordered commensurate 
surface approach a straight line, and the incommensu- 
rate results are between the other two cases. 

The results in Fig. [5] can be understood from a very 
simple mean-field model. The total potential for each 
substrate atom is approximated by an effective harmonic 
spring fc c ff binding it to the substrate plus the potential 
from the rigid atoms. Both are assumed to depend only 
on z, with repulsion from the rigid atoms when z > 
and the equilibrium position relative to the substrate at 
z = —zq. Then 



U(z) = -k eS {z + zof + U w {z) 



(6) 



where U w (z) is the wall potential. For the commensurate 
case, all atoms see the same environment. In the nonad- 
hesive case, U w (z) is just the LJ interaction with the 
rigid atom above, and in the adhesive case it can be cal- 
culated by including the four next-nearest neighbors. In 
principle, a more complicated interaction should be used 
for amorphous and incommensurate simulations because 
of variations in lateral position, but we will see that the 
exact form of U w is not too important. 

For a sufficiently large system we can replace the dis- 
tribution over atoms at a given instant by the equilibrium 
Boltzmann distribution for a single atom. The probabil- 
ity for an atom to be at height z is then 



p{z) 



1 



-U(z)/k B T 



(7) 



9 




0.0 0.2 0.4 0.6 0.8 1.0 



N(f t )/N(0) 

FIG. 5: (Color online) The fraction F(f t )/F(0) of load carried 
by atoms with forces bigger than a threshold ft as a function 
of the fraction of contacting atoms with these large forces 
N(ft)/N(0) for (a) commensurate, (b) incommensurate and 
(c) amorphous surfaces. The threshold f t increases as one 
moves down and to the left along each curve. Circles and 
pluses are for nonadhesive surfaces at ksT/e — 0.175 with 
loads of L/AqE* — 2 x 10" 5 and 0.007, respectively. Triangles 
are for k B T/e = 10~ 4 and load L/A E* = 0.007. Squares 
and diamonds are for adhesive surfaces at k B T/e = 0.175 and 
10" 4 with L/AqE* = 0.007. The dashed lines in (a)-(c) are 
solutions of Eqns. [9] and [10] for adhesive interactions and 
k B T/t = 0.175. The solid line in (a) is a similar solution for 
nonadhesive surfaces. The solid lines in (b) and (c) show the 
simple analytic approximation y — x(l — In a), which is very 
close to the full solution (the solid line in (a)). 

with 

/+oo 
e -U(z)/k B T dz _ (8) 
-oo 

Each threshold force ft corresponds to a threshold height 
Zt, so one can write 

N(f t ) = N / p(z)dz, (9) 

and 

/■oo 

F(ft) = N / f(z)p(z)dz, (10) 



where No is the total number of substrate surface atoms. 
These integrals are readily performed numerically. Re- 
sults for all systems can be fit by varying zq at a fixed 
k c ff — 86e/c 2 . This value of k e s is slightly smaller than 
the estimate of 2k obtained below Eq. [3] by considering 
fixed nearest-neighbors, which reflects the extra compli- 
ance associated with the motion of these neighbors. 

Note that all surfaces follow nearly the same behavior 
at high temperatures, and amorphous results are simi- 
lar even at the extremely low temperature of 10~ 4 e/fcB- 
This surprisingly universal asymptotic behavior can be 
understood from an even simpler approximation for small 
fractional contact area. In this limit, P{z) decays rapidly 
with z. Most contacting atoms are near z = and the 
integrals in Eqs. [9] and [10] are dominated by the con- 
tribution near Zt- For small z, the potential from the 
rigid wall can be approximated by a quadratic function 
U w {z) « k'z 2 /2 and f(z) « k'z. In this limit, U w is 
also much smaller than the potential from the substrate, 
and the total potential can be approximated by a linear 
function 

U(z) w kcsZy/2 + kcffZoz = k cS ZQ/2 + z Q k cS f/k' . (11) 

Eq. [9] can then be reduced to N(f t ) = 
c J*?° df exp(— a/) = (c/a) exp(-a/ t ), where 
a = zo^off / 'k' 'k^T . One also has F(f t ) — 
cf™dffexp(-af) = -dN(f t )/da = (l/o + ft)N(F t ). 
Finally, the fractions become N(f t )/N(0) = exp(-aft) 
and F(f t )/F(0) = (1 + af t ) exp(-a/ t ). The resulting 
solid curves in Figs. [5] (b) and (c) have the form 
y = x(l — In a;) and are independent of a. They differ 
only slightly from the numerical integration of F and N 
at high T shown by dashed lines in each panel and the 
solid line in (a). 

The key assumptions in the above approximation are 
that the probability P(z) decays exponentially for z > 
and that the force / = k'z. Fig. [6] shows that the distri- 
bution of instantaneous local forces is exponential for a 
wide range of circumstances. The dependence of the dis- 
tribution on k! , load and temperature has been removed 
by normalizing by the mean force (/) obtained by av- 
eraging over all atoms with an instantaneous repulsion. 
The probability P(f/ (/)) of a given multiple of this mean 
// (/) is then calculated by averaging over all atoms and 
times. 

As expected from Fig. [5] results for all high temper- 
atures follow an exponential form over several decades 
in P. There are deviations at large / because of anhar- 
monicities in the wall potential and variations in atomic 
separation, but these are generally confined to regions 
of very low probability that contribute little to aver- 
age quantities. The simple linear approximation to U{z) 
provides a surprisingly good description over the most 
important force range. Similar results are obtained at 
k-oT/e = 0.07 or about 10% of the melting temperature. 
Thus this exponential behavior should be present at room 
temperature for nearly any solid. 
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FIG. 6: (Color online) Probability density P(f /{/)) as a func- 
tion of the force on an atom normalized by the average force 
//{/}. Data are for three flat upper surfaces: (a) commen- 
surate; (b) incommensurate; (c) amorphous. For nonadhe- 
sive contacts, there are three data sets for k&T/e = 0.175 at 
L/A E* = 2 x 1(T 5 (O), 5.5 x 10 -4 (+), and 0.007 (A), and 
one data set for k B T/e = 10~ 4 at L/A E* = 0.007 (□). For 
adhesive contacts, there is one data set for k B Te — 0.175 and 
L/A E* = 0.007 (x). 



Deviations from exponential behavior increase as T de- 
creases further, and as the load increases. In these limits, 
the mean position zq moves towards and even past z = 0. 
The high load, high temperature behavior is most evident 
for the adhesive cases in Fig. [5] where adhesion leads to 
effective loads that are roughly 5 times higher than the 
largest nonadhcsive loads. The effective mean normal 
pressure L/Aq ~ E* /20, would correspond to ~ lOGPa 
for metals and the fractional contact area is over 50% for 
all surfaces and definitions of contact (Fig. [3]) . Even here, 
the distribution is exponential for 2 decades or more. 

The low temperature limit depends on geometry. For 
the commensurate case, all atoms have the same poten- 
tial and probability. For k B T/e = 10~ 4 and L/A E* = 
0.007, zq is negative and atoms are almost always in 
repulsive contact. As a result, / has a nearly Gaus- 
sian distribution about the mean value. Atoms on the 
incommensurate surface sample different lateral separa- 
tions nearly uniformly. The resulting force distribution 
is also nearly constant. The random height distributions 



on the amorphous surface give something closer to an ex- 
ponential form, but cut off at larger forces. The disorder 
produces something like an effective temperature in this 
case. 

It is interesting to note that exponential distributions 
of forces have also been observed in previous zero temper- 
ature studies of surfaces that are even rougher than our 
amorphous surfaces. Finite-element calculations of self- 
affine fractal surfaces with roughness down to the mesh 
size found a universal exponential decay Q, although 
smoothing surfaces on small scales suppressed the distri- 
bution at low and high pressures [32, 39] . Atomistic stud- 
ies of fractal two dimensional surfaces with lengths up to 
8192 atoms also found an exponential de cay 1251 . Atom- 
istic simulations of 3D rough surfaces 0)0)122, 23] were 
fit instead to a function predicted by Persson's continuum 
theory [1,0] : P(x) = 7rxexp(— irx 2 /4)/2. However this 
analytic form is inconsistent with their observation of a 
monotonic decrease in P with //(/). Their plots of P 
do not appear to be simply exponential either, but were 
only plotted for fractional contact areas of more than 
10%. Their results are discussed further below. 

Another recent study examined contact of amorphous 
surfaces at a temperature of about 8% of the melting tem- 
perature [bl . fl5| . Their amorphous surface was rougher 
than ours, and they found the force distribution was sim- 
ilar to Persson's predicted form. However, the distribu- 
tion decayed less rapidly than Persson's at large forces 
and only the forces from atoms that were in contact 
more than 30% of the time were included in the distribu- 
tion. This choice was motivated by forcing the number 
of atoms contributing to the distribution to equal the 
mean number in contact at any instant [bl . [l5[ . As we 
show next, those atoms that spend least time in contact 
support the smallest force and removing them from the 
average will suppress P at small forces. It would be in- 
teresting to know whether the full force distribution for 
all atoms in Mo et al.'s simulations had an exponential 
form. 

The quantities f c = F(0)/N o and p c = N(0)/N o rep- 
resent the mean force per atom and fraction of time in 
contact for a given separation zq from the wall. For the 
commensurate case, all atoms will have the same separa- 
tion and same average properties. As the load increases, 
zq will decrease, and both f c and p c will increase. Atoms 
on the amorphous and incommensurate surfaces have dif- 
ferent environments, and thus sample different zq at the 
same load. If they act independently, then zq is the only 
relevant parameter and results for f c vs. p c should fall 
on a universal curve at each temperature. 

Fig. [7£a) shows this collapse for nonadhesive surfaces. 
For all systems, results for a range of loads are combined 
to sample the full range of p c . The results are averaged 
over atoms within narrow ranges of p c and the fluctua- 
tions are comparable to the symbol size. The data are in 
excellent agreement with the simple harmonic mean-field 
model, whose predictions for k^T '/e = 0.175 and 0.07 are 
shown by lines. Note that the displacements associated 
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with the largest forces are quite large and are close to the 
threshold for plastic deformation. 
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FIG. 7: (Color online) (a) and (b) Time-average contact force 
on an atom f c = F(0)/Nq vs. the fraction of time it is 
in contact p c = N(0)/Nq. (a) nonadhesive contacts of var- 
ious upper surfaces: commensurate spherical tip (A), incom- 
mensurate spherical tip (□), amorphous sphetical tip (O)j 
stepped spherical tip (v)> commensurate flat surface (+), in- 
commensurate flat surface (x), amorphous flat surface ((>)■ 
(b) adhesive contact of an amorphous flat surface. The tem- 
perature is T — 0.175e/A;B except for the indicated set of data 
at T = 0.07e/A:B. The solid and dashed lines are fits using 
Eqs. 1101 and[9l with the wall potential and fc e ff used in earlier 
figures, (c) Fraction of atoms in contact more than a fraction 
p c of the time as a function of p c for nonadhesive amorphous 
surfaces at loads L/A E* = 2.6 x 1(T 6 , 2 x 1CT 5 , 1.6 x 1CT 4 , 
5.5 x 10" 4 , 2.6 x 10" 3 , and 7 x 10~ 3 from left to right. 

Fig. [7£b) shows similar results for adhesive surfaces 
pfjj . The agreement with simple theory is still quite 
good, but the fluctuations between atoms (errorbars) are 
greater. We find that the tails in the adhesive potential 
from nearby rigid atoms lead to different shapes of U w for 
substrate atoms in different environments, and there is 
also some coupling of lateral and normal displacements. 
While the distribution of force in Fig. [5] is not affected 
by these changes in shape, the total forces and probabili- 
ties are. Note that the distributions for nonadhesive and 
adhesive contacts in Fig. [7]are fairly similar even though 
the LJ interaction is half as strong in the adhesive case. 
The reason is that interactions from multiple atoms are 
usually important for the adhesive case and they roughly 
double the effective strength of the potential. This simple 



approximation is used in the fit line. 

While the relation between f c and p c is independent of 
load, the distribution of values for different atoms along 
a surface changes with load. For the commensurate case, 
all atoms have the same p c and f c in long time aver- 
ages. For amorphous and incommensurate surfaces the 
quenched disorder leads to a distribution of values cor- 
responding to different zq. Fig. [7tc) illustrates how the 
range of p c changes with load for the nonadhesive amor- 
phous surface. The fraction of contacting atoms n(p c ) 
that are in contact more than a fraction p c of the time 
is plotted against p c for different loads. As the load in- 
creases, the maximum values of p c (and thus f c ) increase. 
In all cases, a large fraction of the contacting atoms spend 
a very small fraction (< 1%) of their time in contact. 
This fraction of weak contacts drops from more than 90% 
for the lowest load to about 30% for the highest load. At 
larger p c there is roughly linear drop in n(p c ) with the 
logarithm of p c . 

The large number of weak bonds changes the distri- 
bution of forces from the instantaneous results shown in 
Figs. [5]and[6l By analogy with Fig. [5] we define N c (f t ) 
as the number of atoms with time-averaged force f c force 
greater than a threshold ft and F c (f t ) as the force they 
carry. Fig. [Et^a) shows the fraction of load carried by 
the fraction of atoms with the highest forces. The frac- 
tion rises much more rapidly than for the instantaneous 
forces, and depends upon load because the fraction of 
weak contacts drops with load. The fraction of atoms 
carrying 90% of the the load rises from about 15% for 
the lowest load to 30% for the highest load, as compared 
to 60% for the instantaneous force. Knippenberg et al. 
have also found that a very small fraction of atoms carries 
most of the load in their simulations (l8| . 

Fig. |3{b) shows the probability of an atom having a 
given time average force f c as a function of f c . As in 
Fig. [6l the force is normalized by the mean force over 
all contacting atoms (/ c ). However, the curves are not 
independent of load. In each case, there is a large peak at 
low forces corresponding to the weak contacts, and these 
dominate the value of (/). There is a nearly exponential 
distribution of forces among the atoms that carry the 
majority of the load and Fig. EJa) can be fit reasonably 
well by assuming P(f/ (/)) is an exponential plus a delta 
function at zero force. As noted above, Yang and Persson 
studied rougher surfaces and found a P(f/{f)) that was 
monotonically decreasing. Their data are also reasonably 
described by a strong peak at low forces followed by an 
exponential region. 

The above results have direct implications for the re- 
lation between area and load found in earlier sections. 
For the case of commensurate surfaces, all atoms have 
the same value of f c and p c . The fraction of area in 
contact A c /Aq at a given load is just p c and the load is 
L = f c * Nq. The fit lines in Fig. [2ja) come from the 
same fit formula used in Fig. [5] with no free parameters. 
Making the same small area approximation as in Eq. II H 
one can derive an analytic form for small A c /Aq. One 
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FIG. 8: (Color online) (a) Fraction of the repulsive force 
F c (/t)/F c (0) supported by the fraction A/c(/t)/iV c (0) of atoms 
with time averaged force f c larger than a threshold ft (ft de- 
creases from left to right). From left to right the dimension- 
less loads are L/A E* = 6.9 x 1(T 5 , 1.6 x 1(T 4 , 3.2 x 10~ 4 , 

5.5 x 10" 4 , 1.3 x 10" 3 , 2.6 x 10" 3 , 4.4 x 10" 3 and 7 x 10~ 3 . 
(b) Probability P(fc/{fc)) of atoms having a time-averaged 
force f c normalized by the mean force on all contacting atoms 
{f c }. The loads are L/A E* = 2 x 1CT 5 (O), 5.5 x 1(T 4 (□), 

2.6 x 10~ 3 (A), and 7 x 10 -3 (v)- 
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Since both A c and L vary extremely rapidly with y, the 
value of y changes little over the relevant range of A c . 
This explains why A c rises only slightly less rapidly than 
linearly with load in Fig. [5] at small A c /Aq. The main 



temperature dependence comes from Sz ms which rises as 
y/r, explaining the scaling of L with \/T in Fig. [2ja). 

For incommensurate and amorphous surfaces the sit- 
uation is complicated, because the distribution of zq is 
not known. For sufficiently high temperatures, thermal 
fluctuations are more important than geometrical fluctu- 
ations and A c rises roughly linearly with load at small 
loads. However, some fraction of the substrate atoms are 
so far from wall atoms that they never contact during the 
simulation. This can be incorporated by reducing Nq to 
the number of substrate atoms that can contact. The fits 
in Fig. [2] reduce the number of atoms to 84% and 57% of 
No for incommensurate and amorphous surfaces, respec- 
tively. These fractions are consistent with the number 
of contacting atoms in the limits of large load and time 
interval and provide an excellent fit. The correspond- 
ing predictions for k-^T / e = 10 , would be substantially 
above the actual data. The reason is that many fewer 
atoms are close enough to contact, particularly at low 
load. 



IV. CONTACT BETWEEN A SPHERICAL TIP 
AND A FLAT SUBSTRATE 

The case of a spherical tip can also be analyzed with 
the theory developed in the previous section. The main 
difference is that the curvature of the tip leads to addi- 
tional variations in the separation zq between substrate 
and wall atoms. These variations increase in magnitude 
and importance as the radius of the tip decreases. We 
will consider a relatively small radius, R = 100c. This 
is comparable to the radius of AFM tips, leading to con- 
tacts with relatively small numbers of atoms and rapid 
changes in surface separation. 

Fig. [3 a) includes results for the mean force and frac- 
tion of time in contact from simulations with nonadhesive 
spherical tips. Results for all geometries collapse on the 
same universal curve obtained for flat surfaces. A similar 
collapse with data in Fig. [7£b) is found for adhesive tips. 
Tests of the distribution of instantaneous forces confirm 
that the same independent atom model describes sub- 
strate atoms under all tips, although the small number 
of atoms in the contacts means that results from several 
longer simulations are required to get similar statistical 
accuracy. The main new features of the results presented 
below come from changes in zo due to surface curvature. 

The continuum limit of nonadhesive contact by a 
spherical tip is described by Hertz theory It pre- 
dicts that the region of contact is a circle whose radius a 
scales as 



a 
R 



ZL 



AE*R 2 



1/3 



(15) 



The contact area is then ira 2 and is thus proportional to 
L 2 / 3 , rather than rising linearly with load. The pressure 
in the contact depends on the radius r from the center 
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and is given by: 

p(r)/E* = - (3L/AE*R 2 ) 1/3 ^l-r 2 /c 



(16) 



Note that this dimensionless pressure and the radius are 
proportional to the same dimensionless measure of load 
[Lj E* R 2 ) 1 / 3 . The characteristic force per atom f a = 
p(0)A a also grows as the cube root of load, while it is 
linear in load for flat surfaces. This makes it hard for tip 
simulations to span a wide range of forces, and nanometer 
scale tips tend to produce forces at the large end of Fig. 

m 

Luan and Robbins considered the same spherical tips 
used here in the zero temperature limit (l2i . fl3j . They 
found systematic deviations between the time and an- 
gle averaged pressure p(r) and the Hertz prediction (Eq. 
I16|) . The finite range of interactions and the atomic scale 
roughness smeared the pressure over a larger area, and p 
remained nonzero to radii that were twice the predicted 
a at low loads. Deviations were much larger for stepped 
tips, where the area decreased in discrete jumps as new 
terraces made contact. 

Figure [5] shows the load dependence of the area A ob- 
tained by counting contacting atoms over different time 
intervals as a function of (L/E*R 2 ) 1 / 3 . The Hertz pre- 
diction is also shown by a dashed line. This prediction 
is the same for all tips since their surfaces differ by less 
than a molecular diameter. However, they show rather 
different behavior, particularly in the case of the stepped 
tip. 

As for the flat surface results in the previous section, 
the low temperature behavior is relatively insensitive to 
the averaging time and the results are most strongly in- 
fluenced by atomic geometry. The contact area for the 
stepped tip is equal to that of the bottom terrace at low 
loads and jumps up at high loads when the second ter- 
race contacts. The other tips all show the same scaling as 
the Hertz prediction at low loads, but with different pref- 
actors. As shown by Luan and Robbins, the contact is 
spread over a larger radius than predicted by Hertz the- 
ory [HI [H| • All atoms of the commensurate surface that 
lie within this region make contact, leading to an area 
about twice as large as the Hertz prediction in Fig. [§] 
(a). The quenched variation in substrate-tip separations 
for incommensurate and amorphous surfaces reduces the 
fraction of atoms that contact. This almost exactly com- 
pensates for the increased radius of contact for the in- 
commensurate surface, while the number of contacting 
atoms is roughly half the Hertz prediction for the amor- 
phous surface. The ratio of MD results to Hertz predic- 
tions depends sensitively on the ratio of lattice constants 
in the incommensurate case (25j . and to the density in 
the amorphous case, but is typically within a factor of 
two. 

The time interval is much more important for high 
temperatures. The area corresponding to atoms in con- 
tact at any instant is substantially reduced from Hertz 
theory in all cases. The rise with load is also more 
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FIG. 9: (Color online) Contact area A c vs load L for a spher- 
ical tip with different geometries: (a) commensurate; (b) in- 
commensurate; (c) amorphous; (d) stepped. Open and filled 
symbols are for T = 0.175e/fcB and 10 _4 e/fcB, respectively. 
The contact area is measured by counting the number of 
atoms in the top layer of the substrate that interact with the 
opposite surface at any instant At = (O) or during time 
intervals At — 0.5r (□) or 500r (A). The dashed lines repre- 
sent the Hertz prediction and are the same in all panels. Solid 
lines represent the fits for each tip from the simple harmonic 
mean-field theory developed in Sec. IIII Cl with ksT/e = 0.175 
and No set equal to the number that contact at k^T/e = 1CP 4 . 



rapid than the Hertz prediction, showing something like 
the more linear area-load relation found for flat surfaces. 
Note that the range of dimensionless pressures is much 
higher for tips because of the small radius. The Hertz 
prediction for our geometry and loads gives p/E* be- 
tween 0.01 to 0.07, while for fiat surfaces p/E* — L/A E* 
ranged from less than 10~ 5 to 0.01 in Fig. O 

The solid lines in Fig. (Hla)-Ellc) show the prediction of 
the harmonic mean-field theory for the commensurate, 
incommensurate, and amorphous tips, taking the corre- 
sponding low temperature area as Aq. The results are in 
excellent agreement with the simulations and imply that 
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the contact area of tips should scale linearly with load at 
low p/E*. It is not possible to test this scaling to lower 
loads in our simulations since there are only a handful of 
atoms touching at the two lowest loads. Accessing the 
lowest p/E* studied for flat surfaces would require in- 
creasing R into the micrometer range. This size scale is 
relevant for nanoindcntcrs and is consistent with common 
estimates of the size of surface asperities in macroscopic 
contacts [EEl- 

As the averaging time interval At grows, the contact 
area rises above the Hertz prediction. The increase is 
more than an order of magnitude at low loads. The re- 
sults also become much less sensitive to tip geometry, be- 
cause thermal fluctuations overcome the quenched vari- 
ation in height. The magnitude of the increase in area 
can be estimated from the tip geometry and our esti- 
mate of height fluctuations in Eq. [3] Near the edge of 
a contact of radius a, the height varies with change in 
radius dr as h(r) « h(a) + adr/R. Thus a fluctuation 
in height by dz can allow contact out to a radius that 
is larger by dr w Rdz/a and an area that is larger by 
dA = lixadr = 2irRdz. If we assume that dz is about 
three times the standard deviation Sz rms over the course 
of the simulation, then dA ~ 70a 2 for R = 100a and 
Sz = 0.04(T. The observed changes are comparable to 
these simple estimates. The fractional change in area 
scales as 5z ms / R and would become smaller for microm- 
eter and larger tips. 

The time averaged pressure is small in the region where 
contact is only made possible by thermal fluctuations. 
One consequence is that measures of the contact radius 
based on the second moment of the pressure distribution 
remain closer to the value predicted by Hertz and mea- 
sured at low temperature [26(. Another is that a small 
fraction of the atoms carry a very large fraction of the 
load. This observation is similar to the result for time av- 
eraged forces between flat surfaces at low temperatures, 
but is due to the large ring at the perimeter of the contact 
that is within Sz of the mean substrate height. 

Figure |TD] shows corresponding results for adhesive 
tips. Results for all but the stepped tip look very similar 
to the flat amorphous case. The area rises sublinearly 
near the pulloff force, and the curves become more lin- 
ear at large loads. The time interval has little effect at 
low temperatures because height fluctuations are smaller 
than geometrical features. The instantaneous contact 
area is smaller at high temperatures, but A c grows with 
At as thermal fluctuations bring more atoms into con- 
tact. 

What is intriguing is that while the results for adhesive 
tips look similar to those for flat amorphous solids, they 
are also qualitatively similar to continuum predictions 
for sphere-on- flat contact. The form of the continuum 
theory depends on the work of adhesion per unit area w 
and the range of interactions ho. The key dimensionless 
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FIG. 10: (Color online) Contact area A c vs load L for a spher- 
ical tip with different geometries: (a) commensurate, (b) in- 
commensurate, (c) amorphous, and (d) stepped. Open and 
filled symbols are for T = 0.175e/fc B and 10" 4 e/fc B , respec- 
tively. The contact area is measured by counting the number 
of atoms in the top layer of the substrate that ever repel the 
opposite surface at any instant At = (O) or during time 
intervals At = 0.5r (□) or 500r (A). Also shown is a fit based 
on the M-D theory to the data for T = 10~ 4 e/fc B . 



parameter scales as (42H4 



/ 9Rw 2 



\2nE* 2 hl 



1/3 



(17) 



The limit of strong, short range interactions (A > 5) is de- 
scribed by Johnson-Kendall-Roberts (JKR) theory, while 
the opposite limit (A < 0.1) is described by Derjaguin, 
Mullcr and Toporov (DMT) theory. For typical scanning 
probe microscope tips and our simulations, A ~ 0.1 to 
1 lies between JKR and DMT limits. Maugis-Dugdale 
(M-D) theory provides a description of this intermediate 
region, and approximate analytic expressions for this the- 
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ory have been developed to simplify fits to M-D theory 

put 

Except for the stepped tip, the data shown in Fig. [TO] 
can be fit reasonably well to M-D theory if the pulloff 
force L p and work of adhesion w are treated as fit pa- 
rameters. Fit parameters are summarized in the Table 
I. While the fits look reasonable, the surfaces often sep- 
arate at forces that are significantly less negative than 
the fit L p . This is particularly evident for FigfTUl where 
the system was not stable at any value lower than the 
final data point. In addition, the fit values of w are all 
larger than the independent measures of w obtained pre- 
viously for the same surfaces at zero temperature (cap- 
tion) [T3| . Presumably this increase in the fit w is needed 
to compensate for the increase in contact area due to the 
atomistic effects discussed above. The fit parameters are 
generally within a factor of two of independently mea- 
sured quantities and may thus be useful in extracting 
approximate values from experiments. However, our re- 
sults reinforce the conclusion reached previously [T2T - fl4j 
that the success of fits to M-D theory does not represent 
a quantitative success of continuum theory, which can 
not be expected to hold at atomic scales. 

The contact areas in Fig. [TOl were obtained by counting 
all atoms that ever felt a repulsion during the given At. 
Fig. [TT] shows the area obtained by counting only atoms 
that have a net repulsive force after averaging over At. 
As in the case of contact between two flat surfaces, this 
definition of A c is insensitive to both averaging time in- 
terval and temperature. It is also harder to describe with 
continuum theory. The area grows smoothly with load, 
but could not be fit to M-D theory. The reason is that 
the power law describing the rate of increase in area is 
significantly different from the value of 2/3 that enters 
continuum theories. For flat incommensurate and amor- 
phous surfaces the area rose smoothly rather than tran- 
sitioning rapidly to full contact. A similar effect seems 
to modify the rise in area with load for the case of spher- 
ical tips. As for flat surfaces the suppression is larger for 
amorphous surfaces than incommensurate surfaces. 



V. DISCUSSIONS AND CONCLUSIONS 

The results presented in this paper reveal the diffi- 
culties in extending the notion of contact to molecular 
scales, even for the simplest geometries. In the case of 
parallel, atomically flat surfaces, continuum theory would 
predict a sharp transition from no contact to full contact. 
The transition would occur at zero load for nonadhesive 
surfaces, and also for adhesive surfaces if contact is iden- 
tified with repulsion. This sharp transition is only seen 
in atomic simulations for the highly unlikely case of two 
identical, commensurate surfaces with atoms aligned on 
top of each other. Even there, it is only found near zero 
temperature or using time averaged forces with adhesive 
interactions. In all other cases, the contact area rises 
very gradually. Indeed, pressures comparable to the ideal 
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FIG. 11: (Color online) Contact area A c based on average 
force (contact means average force being repulsive) vs load L 
for different tip geometries: commensurate (□), incommensu- 
rate (A), and amorphous (O)- Open (filled) symbols are for 
T = 0.175e/fc B (l(T 4 e/M- 



hardness are required to reach full contact. 

Measures of contact area based on instantaneous forces 
between atoms were shown to be extremely sensitive 
to thermal roughness. One common approach is to 
identify the contact area A c with the mean number of 
atoms feeling a repulsive force times the area per atom 
[li . flBl liH I251 [3fjj . For both flat surfaces and spherical 
tips, thermal fluctuations lead to a linear rise in A c with 
load when the mean pressure in the contact is small. For 
flat commensurate surfaces, linear scaling extends over 
the entire temperature range studied, corresponding to 
T from 0.01% to 25% of the melting temperature. For 
incommensurate and amorphous surfaces, the linear be- 
havior is suppressed at low temperatures, as the variation 
in local geometry around surface atoms becomes more 
important than thermal fluctuations. At ksT/e = 10~ 4 , 
A c rises sublinearly, but still remains below full contact 
until very high loads. 

A simple mean-field model was developed that quan- 
titatively describes changes in A c . Each substrate atom 
is assumed to be bound to the substrate by a harmonic 
spring k c ff that represents the net effect of bonds to neigh- 
boring substrate atoms. The atom also feels a force from 
the opposing wall U w (z) that depends on its height z. 
The probability of an atom having a given height and 
force can then be calculated from the Boltzmann equa- 
tion. Over a wide range of parameters only the effective 
stiffness k' associated with the second derivative of U w 
near its minimum is important. 

In addition to reproducing the variation of A c with 
load, this simple harmonic mean-field model makes sev- 
eral predictions about local forces that are consistent 
with simulations. One is that the distribution of instanta- 
neous forces / on substrate atoms is exponential, which 
leads to a universal distribution of the fraction of load 
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TABLE I: Parameters from fits of contact area in Fig. [10] to M-D theory. Independently determined values for wa 2 /e at zero 
temperature are 1.05, 0.45 and 0.23 for commensurate, incommensurate and amorphous cases. 

tip parameter T = 0.175e/fc B T = 10" 4 e/fc B 







0.005r 
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commensurate 


L p {e/a) 


592 
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558 
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w(e/a 2 ) 


1.11 


1.21 


1.38 


1.19 


incommensurate 


L p (e/a) 


480 


477 


483 


493 




w(e/a 2 ) 


0.84 


0.93 


1.08 


0.90 


amorphous 


L p (e/a) 


391 


331 


332 


417 




w{e/a 2 ) 


0.65 


0.64 


0.79 


0.70 



born by the atoms with the largest forces. Figs. [5] and [6] 
show that these predictions describe contact between flat 
surfaces with all geometries except at extremely low tem- 
peratures (< T m /10) and high pressures (p/E* > 0.01). 

Another prediction is that for all atomic (commensu- 
rate, incommensurate, amorphous, or stepped) and large 
scale (flat or sphere) geometries, the mean force on an 
atom and the fraction of time in contact have a func- 
tional relationship that depends only on temperature. 
The same relation holds for adhesive surfaces if contact 
is defined as repulsion and the mean repulsive force f c is 
calculated. The fraction of time in contact rises linearly 
with force at low forces and this quantitatively describes 
the variation of A c with load for commensurate surfaces. 
For incommensurate and amorphous surfaces, the predic- 
tion need only be scaled by a constant factor to include 
the fact that some substrate atoms are too far from wall 
atoms to make contact. The linear relation between A c 
and load or p c and f c breaks down when the fraction of 
time in contact approaches unity. This does not occur 
until quite large values of p/E* . For example, p c reaches 
50% at p/E* = 0.01 and 0.02 for T/T m = 0.1 and 0.25, 
respectively. 

The simplicity of the harmonic mean- field model allows 
us to make estimates about the range of validity of the 
above statements that are quite general and should not 
depend on the details of atomic interactions. Equation 
[13] gives the ratio between the dimensionless pressure and 
the fraction of area in contact at any instant. The main 
variation in this ratio comes from <5z rms , which rises as 
the square root of temperature. The Lindemann criterion 
says that the value of Sz rms should be about 10% of the 
nearest- neighbor spacing at the melting temperature . 
This allows us to write 5z rms « O.lay^T /T m . Inserting 
this into Eq. [T3]and multiplying and dividing by k e g we 
find 

A c L 



A CA A E* 



= CA 



_P_ 

E* 



with 



ca 



ln E*A a k e g 
10— rrV 



(18) 



(19) 



crkes k' 

The values of E* and k e g are both determined by the 
interactions between neighbors, and the ratio E* A a /ak e g 



will generally be near unity. For the potential used here 
it is 0.92. The value of k'/k e g should also be close to 
unity if interactions across the interface are comparable 
to internal interactions. The value of y is directly related 
to A c /Aq via Eq. [T^l but varies extremely slowly. For 
example, it drops from 3.1 to 1.4 as A c /Aq rises from 
10 -3 to 10" 1 . One can thus conclude that 



ca 



(20) 



for any potential. Note that this expression differs from 
the simple estimate based on thermal collisions in Eq. [5] 
by a factor of order unity. 

Molecular dynamics is only accurate for temperatures 
high enough compared to the Debye temperature that 
quantum effects can be ignored. This typically corre- 
sponds to T/T m > 0.05 and applies to most materials at 
room temperature. As T increases from 0.05T m to T m 
the pressure required for atoms to be in contact more 
than 50% of the time rises from p/E* ~ 0.01 to 0.05. 
Plastic deformation sets in when p exceeds the hardness, 
H. Typical values of H/E* are of order 10~ 4 to 10~ 2 
for macroscopic crystals. While larger values may be ob- 
served in defect free systems and at nanometer scales, the 
theoretical limit is only H/E* ^0.1. Thus one can gen- 
erally expect that atoms will spend a significant fraction 
of their time out of repulsive contact at temperatures and 
pressures of experimental interest. 

Our simulations reduced the number of free parame- 
ters and the size of the simulation by keeping one surface 
rigid. If both surfaces are compliant, there are two offset- 
ting changes. One is that there are thermal fluctuations 
by atoms on both sides that increase thermal roughness. 
This can be modeled roughly by replacing fc^ 1 by the 
sum of the inverse values for the two surfaces. Contin- 
uum contact mechanics says that the inverse of the ef- 
fective modulus is related in the same way to the moduli 
of the two solids. Since the same interactions determine 
both E* and k e g, the two changes should nearly cancel. 
Simulations for a few examples confirmed this. 

The above discussion neglects quenched geometric dis- 
order due to differences in separation between substrate 
atoms and atoms on the opposing surface. This reduces 
the number of atoms on amorphous and incommensu- 
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rate surfaces that can be brought into contact by thermal 
fluctuations. At T/T m = 0.1 and 0.25 we found that the 
mean-field theory still worked if one reduced the available 
number of contacts by a constant fraction. For example, 
the fit lines in Fig. 5 used 84% for the incommensurate 
surface and 57% for the amorphous case. 

At the much lower temperature of 10~ 4 e//cB, Sz Tma ~ 
0.001<7 is much smaller than the atomic scale roughness 
on incommensurate and amorphous surfaces. Particu- 
larly in the case of amorphous surfaces, those atoms that 
contact at a given load stay in contact most of the time. 
This is clearly seen in Fig. 2, where for all loads A c (At) 
increases by less than a factor of 2 as At/r increases 
from to 500. One can also calculate the mean force 
on atoms from the measured A c and L. At the lowest 
load in Fig. 2c, f c ~ 0.2e/er, which is larger than the 
force k'5z rms ~ O.le/er where the linear relation between 
force and p c breaks down. Thus, in contrast to the com- 
mensurate case, the rise in A c reflects an increase in the 
number of atoms that are close enough to contact rather 
than an increase in the fraction of time in contact. As 
load increases, contacting atoms are pushed down rela- 
tive to other atoms, allowing new atoms to contact the 
opposing surface. Note that the rise in A c with load is 
nonlinear for both incommensurate and amorphous sur- 
faces. Continuum theory predicts a linear relation for 
rough surfaces, but it is not clear it can be applied at 
these scales or that the surfaces have an appropriate dis- 
tribution of roughness. 

Experimental surfaces are generally much rougher than 
those considered here, with bumps on top of bumps at 
all scales [l|, At low loads, contact occurs only where 
two asperities from opposing surfaces overlap. In contin- 
uum theories, the linear area- load relation comes from 
the growth in the number of such contacts with load. 
The distribution of contact sizes and forces remains un- 
changed. 

Simulations of the spherical tip geometry give in- 
sight into the behavior of asperity contacts. As for flat 
surfaces, at typical experimental temperatures, thermal 
roughness leads to a linear relation between area and 
pressure until p/E* is 0.01 to 0.05. This pressure is gen- 
erally large enough to produce plastic deformation under 
micrometer and larger asperities, and comparable to the 
hardness of nanoasperities. In the cases of nanoasperi- 
ties, these and previous simulations [III, [U, [26| show that 
the contact area is not accurately described by continuum 
theory. Thus it is not clear that A c will follow contin- 
uum predictions [U, [E] for rough surfaces at experimental 



temperatures. 

Since continuum theory ignores thermal fluctuations, 
it is natural to work with definitions of contact area based 
on time-averaged rather than instantaneous forces. For 
adhesive surfaces, the area where time-averaged forces 
are repulsive is nearly independent of averaging time 
and is insensitive to temperature. While this seems the 
least ambiguous definition of contact, it only shows a 
sharp transition from no contact to full contact for flat 
commensurate surfaces. For flat incommensurate and 
amorphous surfaces, the transition is spread over a range 
Ap/E* ~ 0.03 that is comparable to the pressures needed 
to produce plasticity. 

The time averaged force gives much less satisfactory re- 
sults for nonadhesive surfaces. Since the force is always 
repulsive, any contact leads to a positive time-averaged 
force. The contact area grows monotonically with the 
averaging time interval At, leading to substantial ambi- 
guity in A c for both flat and spherical surfaces. Aver- 
aging positions instead of forces reduces the sensitivity 
to rare events. However, one finds that the mean sep- 
aration is beyond the interaction range for most loads, 
implying that there is no contact even though there is a 
large repulsive force. 

Averaging positions and forces gives similar results for 
adhesive surfaces if contact is based on repulsive interac- 
tions between atoms. However we have found that an- 
other definition gives a sharper transition from partial 
to full contact for incommensurate surfaces. Instead of 
basing contact on forces between atoms, it is based on 
the separation between surfaces. A Delaunay triangula- 
tion of each surface is performed and triangles on oppos- 
ing surfaces are said to contact if they are separated by 
less than the distance where atomic interactions become 
repulsive. We are currently exploring whether associat- 
ing contact with surface separations instead of atomic 
forces improves the comparison of simulations and con- 
tinuum calculations of rough surfaces with multiaspcrity 
contacts. 
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